#ifndef PCPS_PROPERTIES_HEADER
#define PCPS_PROPERTIES_HEADER

#include <src/pcps.h>
#include <src/input.h>
#include <src/hamiltonian.h>
#include <src/wavefunction.h>
#include <src/accumulator_properties.h>
#include <src/stochastic_optimization.h>

namespace pcps {

  //--------------------------------------------------------------------------------
  // pcps::properties -- Computes correlation functions and writes them to a file.
  //--------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF>
  inline void properties(const pcps::Input & userinp,
                         const pcps::Hamiltonian<SCALAR> & ham,
                         pcps::Wavefunction<SCALAR> & wfn,
                         long int sample_length, // input is the TOTAL number of samples from all threads
                         const std::string & filename) {

    // get MPI info
    const MPI::Comm & comm = MPI::COMM_WORLD;
    const int nproc = comm.Get_size();
    const int myrank = comm.Get_rank();

    // divide the sample length among the different threads
    sample_length = 1 + sample_length / ( nproc * userinp.nthreads() );

    // define a short name for a shared pointer to a process accumulator
    typedef boost::shared_ptr<pcps::ProcessAccum<SCALAR, PART, REF> > ProcAccPtr;

    // initialize a vector of accumulators and place the properties accumulator in it
    std::vector<ProcAccPtr> accumulators;
    { ProcAccPtr temp( new pcps::PropertiesProcessAccum<SCALAR, PART, REF>(userinp, wfn, sample_length) ); accumulators.push_back(temp); }

    // compute the properties via stochastic sampling
    pcps::stochastic_iterations<SCALAR, PART, REF>(userinp, ham, wfn, accumulators, 1, 2, 2, sample_length, userinp.burn_in_length());

    // write the computed correlation functions to a file
    if (myrank == 0) ( dynamic_cast<pcps::PropertiesProcessAccum<SCALAR, PART, REF> *>(accumulators.at(0).get()) )->write_to_file(filename);

  }

} // end namespace pcps

#endif
